Phase diagram and thermodynamics of the three-dimensional Bose-Hubbard model 
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We report results of quantum Monte Carlo simulations of the Bose-Hubbard model in three 
dimensions. Critical parameters for the superfluid-to-Mott-insulator transition are determined with 
significantly higher accuracy than it has been done in the past. In particular, the position of the 
critical point at filling factor n — 1 is found to be at (U/i) c ~ 29.34(2), and the insulating gap A is 
measured with accuracy of a few percent of the hopping amplitude t. We obtain the effective mass 
of particle and hole excitations in the insulating state — with explicit demonstration of the emerging 
particle-hole symmetry and relativistic dispersion law at the transition tip — along with the sound 
velocity in the strongly correlated superfluid phase. These parameters are the necessary ingredients 
to perform analytic estimates of the low temperature (T <C A) thermodynamics in macroscopic 
samples. We present accurate thermodynamic curves, including these for specific heat and entropy, 
for typical insulating (U/t = 40) and superfluid (t/U = 0.0385) phases. Our data can serve as a 
basis for accurate experimental thermometry, and a guide for appropriate initial conditions if one 
attempts to use interacting bosons in quantum information processing. 
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I. INTRODUCTION 

In the past decade, strongly correlated lattice quan- 
tum systems have been attracting a lot of interest and 
effort. Remarkably, simple yet nontrivial models which 
contain most of the important many-body physics and 
known in the theory community for many years can be 
now realized and studied experimentally. For the first 
time theoretical predictions and experimental data for 
strongly correlated states can be directly tested against 
each other in the ideal setup when all model ingredients 
are known and controlled. 

Experimentally, lattice systems are realized by trap- 
ping atoms in an optical lattice, a periodic array of po- 
tential wells resulting from the dipole coupling of the 
atoms to the electric field of the standing electromag- 
netic wave produced by a laser. Optical lattices are a 
very powerful and versatile tool. By changing the laser 
parameters and configuration, the properties and geom- 
etry of the optical lattice can be finely tuned Ul- 
timately, this results in the possibility of controlling the 
Hamiltonian parameters and exploring various regimes of 
interest. In particular, ultra-cold Bose atoms trapped in 
an optical lattice are an experimental realization of the 
Bose-Hubbard model. The model has been studied in 
the seminalpaper by Fisher, Weichman, Grinstein, and 
Fisher, Ref.y, and its physical realization with ultra-cold 
atoms trapped in an optical lattice has been envisioned 
m Ref. 31. Few years later, the Bose-Hubbard system 
was produced in the laboratory EL S ince then, the field 
remains very active p| H 0, H, l9l. Ilfjj] , not only because 
theoretical predictions and experimental techniques still 
have to be substantially improved to claim the quanti- 
tative agreement, but also because of the new physical 
applications. 

At zero temperature, a system of bosons with com- 



mensurate filling factor undergoes a superfluid-to-Mott 
insulator (SF-MI) quantum phase transition. The ground 
state of MI can be used in quantum information process- 
ing to initialize a large set of qubits (the main remaining 
challenge is in addressing single atoms to build quantum 
gates, see Ref. [l[ and references therein). Atomic sys- 
tems in optical lattices have the advantage of being well 
isolated from the environment. This results in a relatively 
long decoherence time of the order of seconds [l[ and 
therefore the possibility of building long-lived entangled 
many body states. These properties make MI ground- 
states good candidates for building blocks of a quantum 
computer. Another possible application is in interfero- 
metric measurements . It has been argued [HI, [l3|, [3] 
that using the superfluid-to-Mott-insulator phase transi- 
tion to entangle and disentangle atomic Bose-Einstein 
condensate one can go beyond the Heisenberg-limited in- 
terferometry. 

A system of bosons with short-range repulsive pair 
interaction trapped in an optical lattice is described by 
the Bose-Hubbard Hamiltonian: 

H = -t h \ b j ■ + \ ruirii - 1) - ^2 ^ > C 1 ) 

<ij> i i 

where b\ and bi are the bosonic creation and annihilation 
operators on the site i, t is the hopping matrix element, 
U is the on-site repulsion and fa = jj,— V(i) is the sum of 
the chemical potential /i and the confining potential V(i). 
In what follows, we consider bosons in the simple cubic 
lattice. At zero temperature and integer filling factor, the 
competition between kinetic energy and on-site repulsion 
induces the MI-SF transition. When the on-site repulsion 
is dominating, t/U <C 1 the atoms are tightly localized in 
the MI ground state which is well approximated by the 
product of local (on-site) Fock states. The Mott state 
is characterized by zero compressibility originating from 
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an energy gap for particle and hole excitations. When 
the hopping amplitude is increased up to a certain criti- 
cal value (t/U) c , particle derealization becomes energet- 
ically more favorable and the system Bose condenses. In 
the chemical potential vs. hopping matrix element plane 
(energies are scaled by U), the T = phase diagram has 
a characteristic lobe shape 0, see also Fig. [3] below, 
with the MI phase being inside the lobe (there is one 
lobe for each integer filling factor). The most interesting 
region in the phase diagram is the vicinity of the lobe tip, 
(p = fj, c , t = t c ), corresponding to the MI-SF transition 
in the commensurate system. For other values of [i or 
t, the SF-MI criticality is trivial and corresponds to the 
weakly interacting Bose gas at vanishing particle density 
Q . It is straightforwardly described provided the particle 
(hole) effective mass is known. If, however, one crosses 
the MI-SF boundary at constant commensurate density 
(this is equivalent to going through the tip of the lobe 
at a fixed chemical potential) the long- wave action of the 
system becomes relativistic and particle-hole symmetric. 
Now the phase transition is in the four-dimensional U(l) 
universality class [2|]. It is worth emphasizing that here 
we have a unique opportunity of a laboratory realization 
of the non-trivial relativistic vacuum, a sort of a "hy- 
drogen atom" of strongly-interacting relativistic quan- 
tum fields. Approaching the critical point from the MI 
side, one deals with the vacuum that supports massive 
bosonic particles and anti-particles (particles and holes). 
On the other side of the transition, the SF vacuum sup- 
ports massless bosons (phonons) that do not have an 
anti-particle analog. In principle, one can systematically 
study universal multiparticle scattering amplitudes of the 
relativistic quantum field theory in the ultra-cold " super- 
collider" ! 

The present study is focused on the three- 
dimensional (3D) system. To the best of our knowledge, 
previous systematic studies of the 3D case were limited to 
the mean- field (MF) 0] and perturbative methods [l5| . 
In Ref. [HI, the authors utilized the strong-coupling ex- 
pansion to establish boundaries of the phase diagram in 
the {n/U, t/U) plane. This approach, based on the small 
ratio zt/U < 1, where z = 6 is the coordination number 
for the simple cubic lattice, works well only in the MI 
phase in the region far from the tip of the lobe, where 
the insulating gap is larger then hopping, A/zt > 1. 
Close to the critical region, where A/zt ~ 1, the strong- 
coupling expansion is no longer valid. We present the 
results of large-scale Monte Carlo (MC) simulations of 
the model (fTJ) by worm algorithm [la ]. With precise 
data for the single-particle Green function, we are able 
to carefully trace the critical and close-to-critical behav- 
ior of the system, and, in particular, produce an ac- 
curate phase diagram in the region of small insulating 
gaps A<t. Though the corresponding parameter range 
is quite narrow, it is crucial to clearly resolve it to re- 
veal the emerging particle-hole symmetry and relativis- 
tic long-wave physics at the tip of the MI-SF transition. 
We also present data for the effective masses of particle 



and hole excitations inside the insulating phase. Close 
to the MI lobe tip, the data for the dispersion of the el- 
ementary excitations are fitted by the relativistic law, in 
agreement with the theory (this also allows us to extract 
the value of the sound velocity in the critical region). 
In the Mott state, the knowledge of gaps and effective 
masses is sufficient to calculate the partition function in 
the low temperature limit analytically and to make reli- 
able predictions for the system entropy. 

For such applications of the system as quantum in- 
formation processing and intcrferomctry, controlling the 
temperature is of crucial importance. Most applications 
are based on the key property of the good insulating 
state, which is small density fluctuations in the ground 
state. At zero temperature fluctuations are of quan- 
tum nature and can be efficiently controlled externally 
through the t/U ratio. At finite temperature, fluctua- 
tions are enhanced by thermally activated particle-hole 
excitations. Only when the temperature is much smaller 
than the energy gap, the number of excitations is expo- 
nentially small. Up to date, there are no available ex- 
perimental techniques to measure the temperature of a 
strongly interacting system. For weakly interacting sys- 
tems, the temperature can be extracted in a number of 
ways, e.g. from the interference pattern of matter waves 
[ItJ or the condensate fraction observed after the trap 
is released and the gas expands freely [l8[ — these prop- 
erties are directly related to the momentum distribution 
function n(k). For strongly interacting systems, both 
temperature and interaction are responsible for filling the 
higher momentum states, which makes it hard to extract 
temperature using absorption imaging techniques. 

The results presented in this paper can be used 
to perform accurate thermometry. Typically, the ini- 
tial temperature, T^ in \ (before the optical lattice is adi- 
abatically loaded) is known. By entropy matching one 
can easily deduce the final temperature of the MI state, 
T = T^ n \ provided the entropy of the MI phase is 
known. To this end we have calculated the energy, spe- 
cific heat and entropy of the system in several impor- 
tant regimes which include MI and strongly correlated 
SF phases. These data can be used to suggest appropri- 
ate initial conditions which make the Bose- Hubbard sys- 
tem suitable for physical applications, such as the ones 
described above. 

Another interesting question concerns the nature of 
inhomogeneous states in confined systems when the MI 
phase is formed in the trap center. The confining po- 
tential provides a scan in the chemical potential of the 
phase diagram at fixed t/U [3| . As one moves away from 
the trap center the system changes its local state. At 
zero temperature, the density profile of the system can 
be read (up to finite-size effects) from the ground state 
phase diagram. At finite temperature, this is no longer 
possible. In particular, the liquid regions outside of the 
MI lobes could be normal or superfluid, depending on 
temperature. 

So far experimental results have been interpreted 
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by assuming that liquid regions are superfluid, but there 
were no direct measurements or calculations to prove 
that this was the case. [Part of the problem is that 
absorption imaging is sensitive only to n(k), which is 
the Fourier transform of the single-particle density ma- 
trix in the relative coordinate. All parts of the system 
contribute to n(k) and it is hard to discriminate where 
the dominant contribution comes from.] It is almost cer- 
tain that T' fin ) /ri fin ' ) of the strongly correlated system 

is higher then T^ ln ' /T^ m '. Indeed, since the entropy of 
MI at A 3> T is exponentially small, most entropy will 
be concentrated in the liquid regions. At this point we 
notice that the transition temperature in the liquid is 
suppressed relative to the non-interacting Bose gas value 
Tc°^ ~ 3.313 n 2 / 3 /m by both (i) effective mass enhance- 
ment in the optical lattice, m — > l/2ta 2 (here a is the 
lattice constant) , and (ii) strong repulsive interactions in 
the vicinity of the Mott phase, in fact, T c — * at the 
SF-MI boundary. It seems plausible that the MI phase is 
always surrounded by a broad normal liquid (NL) region. 
It may also happen that superfluidity is completely elim- 
inated in the entire sample in the final state. [Strictly 
speaking, at T ^ the MI and NL phases are identical 
in terms of their symmetries and are distinguished only 
quantitatively in the density of particle-hole excitations, 
i.e. in the Hamiltonian (TT|) the finite-temperature MI is 
continuously connected without phase transition to NL, 
see Fig. [TJ For definiteness, we will call NL a normal 
finite-T state which is superfluid at T = for the same 
set of the Hamiltonian parameters.] Fig. [T] shows the 
finite-temperature phase diagram for filling factor n = 1 
(we will discuss how we determine the critical temper- 
ature in Sec. III). The critical temperature goes to zero 
sharply, while approaching the critical point. In the limit 
of U — > the critical temperature is slightly above the 
ideal-gas prediction (T = 5.591* was calculated using the 
tight binding dispersion relation), as expected (see, e.g. 
Ref. [H). 

The paper is organized as follows: in Sec. II we 
present results for the ground state phase diagram and 
effective mass of particle (hole) excitations, at integer 
filling factor n = 1. In Sec. Ill we investigate the ther- 
modynamic properties of the system. We present data 
for energy, specific heat and entropy and calculate the 
final temperature of the uniform and harmonically con- 
fined system in the limit of large gaps. For the case of 
trapped system, we also determine the state of the liq- 
uid at the perimeter of the trap. Brief conclusions are 
presented in Sec. IV. 



II. GROUND STATE PROPERTIES 

This section deals with the results of large-scale 
Monte Carlo simulations for the ground state phase di- 
agram of the Bose- Hubbard system in three dimensions. 
Analytical approaches, e.g. the strong coupling expan- 
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FIG. 1: (Color online). Finite-temperature phase diagram 
at filling factor n = 1. Solid circles are simulation results 
(the line is a guidance for the eye), error bars are plotted. 
T — 5.59U is the critical temperature of the ideal Bose gas 
with the tight binding dispersion relation. At finite, but low 
enough temperature, the MI domain is loosely defined as the 
part of the phase diagram to the right of the gray line. The 
rest of the non-superfluid domain is referred to as normal 
liquid (NL). 

sion, work well in the region where zt/U -C 1 and the 
system is deep in the MI phase. Under these conditions 
the kinetic energy term in the Hamiltonian can be treated 
perturbatively and the unperturbed ground state is a 
product of local Fock states. In Ref. [15| the authors car- 
ried out an expansion, up to the third order in zt/U , for 
the SF-MI boundaries and estimated positions of critical 
points at the tips of the MI lobes (by extrapolating re- 
sults to the infinite expansion order) . Their results agree 
with the mean field solution calculated in Ref. [2| , when 
the latter is expanded up to the third order in zt/U and 
the dimension of space goes to infinity. As already men- 
tioned, this approach starts failing when A ~ zt. Using 
MC techniques we were able to calculate critical param- 
eters and predict the position of the diagram tip with 
much higher accuracy: with the worm algorithm (WA) 
approach the energy gaps can be measured with preci- 
sion of the order of 10 _2 i (2(j. The simulation itself is 
based on the configuration space of the Matsubara Green 
function 

G(i,r) = (T T bl(r)b (0)) , (2) 

which is thus directly available. We utilize the Green 
function to determine dispersion relations for particle and 
hole excitations at small momenta [from the exponential 
decay of G(p, r) with the imaginary time] which directly 
give us the energy gap and effective masses. 

Recall that in the momentum space the Green func- 
tion of a finite size system G(p, r) is different from zero 
only for p = p m = 2n(m x /L x , m y /L y , m z /L z ), where 
L a =x,y,z is the linear system size in direction a (we per- 
formed all simulations in the cubic system with L a = L), 
and m = (rn x ,m v ,m z ) is an integer vector. Using 
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FIG. 2: (Color online). Zero-momentum Green function in 
the Mott phase with the chemical potential fi/U = 0.809, 
slightly below the upper phase boundary. Here we show 
data for the system with N — 10 3 bosons at U/t = 70 and 
T/t = 0.025. In the inset we plot the energy gap A for linear 
system sizes L=10 and L=20. Finite-size errors are within 
the statistical error bars. 



Lehman expansion and extrapolation to the r 
limit one readily finds that 
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The two limits describe single-particle/hole excitations 
in the MI phase. Here Z± and e± are the particle/hole 
spectral weight (or Z-factors) and energy, respectively. 
In the grand canonical ensemble, excitation energies are 
measured relative to the chemical potential. With this 
in mind, calculating the phase diagram of the system is 
rather straightforward. At a fixed number of particles 
N = L 3 and t/U ratio one determines chemical poten- 
tials n± for which the energy gap for creating the parti- 
cle/hole excitation with p = vanishes. The insulating 
gap is given then by A = /i + — For high precision 
simulations of the gap one has to choose the value of \x 
very close to [x± and consider finite, but zero for all prac- 
tical purposes, value of temperature so that the following 
two conditions are satisfied: 



\(Jt -fx±\ < t 



\fj, -(M±\ > T 



(4) 



This is exactly how we proceed. By plotting ln[G(p, r)] 
vs. t we deduce e±(p) from the exponential decay of the 
Green function. A typical example is shown in Fig. [2j We 
use the values of the hopping amplitude t and the lattice 
constant a as units of energy and distance, respectively. 

In Fig. Owe present accurate results for the bound- 
aries of the first Mott insulator lobe. We have done cal- 
culations for systems with linear sizes L — 5, 10, 15, 20. 
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FIG. 3: (Color online). Phase diagram of the first MI-SF 
lobe. Numerical data are shown by open circles. The error 
bars are shown but are barely visible even in the inset. The 
dashed lines and the square represent results of Ref. [l5j] for 
the strong coupling expansion and the extrapolated position 
of the diagram tip, respectively. 



Up to values of t/U ~ 0.031 no size effects were detected 
within the error bars. [Here and throughout the paper 
error bars are of two standard deviations] . In the critical 
region the finite-size effects were eliminated using stan- 
dard scaling techniques (see below). The dashed lines 
are the prediction of Ref. [lj| based on the third-order 
expansion in t/U. It becomes inaccurate quite far from 
the tip when the insulating gap is about ~ 6t. On the 
other hand, the value of the tip position extrapolated 
to the infinite order is right on target, within the error 
bar of order 3t for the chemical potential and on-site 
repulsion. In all simulations (performed at t/T = 40) 
the finite-temperature effects are negligible — the system 
is essentially in its ground state. 

To eliminate finite-size effects in the critical region 
and pinpoint the position of the lobe tip, we employed 
standard scaling techniques based on the universality 
considerations. 

First, let us briefly review the universal proper- 
ties of the insulator-to-superfluid transition (see Ref. 
for more details). There exist two types of transi- 
tions: the "generic" transition, when the phase bound- 
ary is crossed at fixed t/U, and a special transition 
at fixed integer density, when the SF-MI boundary is 
crossed at fixed fi/U. The generic transition is driven by 
the addition/subtraction of a small number of particles, 
and is fully characterized by the physics of the weakly- 
interacting Bose gas formed by the small incommensurate 
density component n — no, where no is the nearest integer 
to n. In particular, if S is the deviation from the generic 
critical point in the chemical potential or t/U ratio then 
\n — no | - <5 and T C (S) - S 2 / 3 in the SF phase. 



5 



The special transition at the tip of the lobe hap- 
pens at fixed integer density. It is driven by delocalizing 
quantum fluctuations which for large values of t/\J en- 
able bosons to overcome the on-site repulsion and hop 
within the lattice. As explained in Ref. 0, the effective 
action for the special transition belongs to the (d + 1)- 
dimensional XY universality class which implies emer- 
gent relativistic invariance (rotational invariance in the 
imaginary-time-space, which is equivalent to the Lorentz 
invariance in real-time- space), and, in particular, an 
emergent particle-hole symmetry. The upper critical di- 
mension for this transition is (d + 1) = 4, so that for 
d > 3 the critical exponents for the order parameter, j3, 
and the correlation length, v, are of mean-field character: 
j3 = v = 1/2 (with logarithmic corrections for d = 3). In 
this study, we were not able to resolve logarithmic renor- 
malizations for realistic 3D systems and proceed below 
with the analysis which assumes mean-field scaling laws. 

Denoting the distance from the critical point as 7 = 
[(t/U) c — t/U], for a system of linear size L one can write 

A( 7 ,£) = rV(£/£) = L-'g^L 2 ), (5) 

where A is the particle-hole excitation gap, £ is the corre- 
lation length, and f(x) and g{x) are the universal scaling 
functions. In the last expression we have used the rela- 
tion £ oc 7" 1 / 2 . At the critical point, the product LA 
does not depend on the system size. Therefore, by plot- 
ting LA as a function of t/U one determines the critical 
point from the intersection of curves referring to differ- 
ent values of L, as shown in Fig. 0J This analysis yields 
(Fig. [5] explains how finite-size effect in the position of 
the crossing point originating from corrections to scaling 
was eliminated) 

(t/U) c = 0.03408(2) (n = 1) . (6) 

Our final results are summarized in Fig. [3] We find 
that the size of the critical region where 4D XY scal- 
ing laws apply is narrow and restricted to small gaps of 
the order of A < t (inside the vertical error bar on the 
strong-coupling expansion result in Fig. [3]). It appears 
that resolving this limit experimentally would be very 
demanding. 

To perform analytic estimates of the MI state en- 
ergy (and entropy) at low temperature T <C A one has 
to know effective masses of particle and hole excitations, 
m±. For example, the particle/hole contributions to en- 
ergy in the grand canonical ensemble are given by the 
sums 

E ± = £e±(k)n e « (^) 3 / dke ± (k)e-^ T , 

where n e is the Bose function and 

e±(k) « ±( M± - n) + k 2 /2m ± (8) 
For large gaps the tight binding approximation 

e ± (k) * ±( M± -/,)+ £ (Q) 

a—x,y,z 




FIG. 4: (Color online). Finite size scaling of the energy gap at 
the tip of the lobe. AL/t vs. t/U for system size L=5 (solid 
squares), L=10 (open circles), L—15 (solid circles), L=20 
(open squares). Lines represent linear fits used to extract 
the critical point. 




FIG. 5: Extrapolation to the thermodynamic limit. We 
show the intersections (triangles) of the curves (L=5, L=10), 
(L=10, L=15), (L=15, L=20), vs. L~l x . The fit (solid line) 
yields {t/U) c = 0.03408(2). 

is a reasonable approximation for all values of k in the 
Brillouin zone. Note that, if one is to use the local den- 
sity approximation (LDA) for the energy/entropy esti- 
mates of trapped systems, then calculations have to be 
performed in the grand canonical ensemble. 

To determine effective masses we computed G(p, r) 
in the insulating state and deduced e±(p) for several low- 
est momenta from the exponential decay of the Green 
function on large time scales. Dispersion laws were then 
fitted by a parabola, with the exception for the diagram 
tip, where the dispersion relation is relativistic. The re- 
sult for m± is shown in Fig. [S] When t/U ^0 one can 
calculate effective masses perturbatively in t/U to get 

tm + = 0.25 -3t/U, tm- = 0.5-12t/U. (10) 

Clearly, our data are converging to the analytical result 
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FIG. 6: (Color online). Effective mass for hole (solid cir- 
cles) and particle (open circles) excitations as a function of 
t/U. The exact results at t/U — are m+ = 0.25/t and 
m_ = 0.5/t. By dashed lines we show the lowest order in t/U 
correction to the effective masses. Close to the critical point 
the two curves overlap, directly demonstrating the emergence 
of the particle-hole symmetry. At t/U = 0.034, the sound 
velocity is c/t = 6.3 ± 0.4. 




FIG. 7: (Color online). Energy per particle at t/U = 0.005, 
unity filling factor and linear system size L — 20. Solid cir- 
cles are prime data (error bars within symbol size), the solid 
line is the analytical prediction from Eq. [7] where, at each 
temperature, the chemical potential has been fixed by impos- 
ing equal number of particle and hole excitations. The inset 
shows the total number of excitations present in the system. 



as t/U — > 0. On approach to the critical point the ef- 
fective mass curves become identical for particles and 
holes indicating that there is an emergent particle-hole 
symmetry at the diagram tip. In agreement with the 
theoretical prediction, the data taken at t/U = 0.034 
are fitted best with the relativistic dispersion relation 
e(p) = CyJ m 2 c 2 + p 2 , where c is the sound velocity and 
the effective mass is defined as m* = A /2c 2 . At this 
value of t/U we have found that c/t = 6.3 ± 0.4 and 
tra» = 0.010 ±0.004. 



III. FINITE TEMPERATURE ANALYSIS 

Controlling the temperature is an important experi- 
mental issue, crucial for many applications of cold atomic 
systems and studies of quantum phase transitions. In 
this section, we discuss thermodynamic properties of the 
Bose-Hubbard model. We present data for energy, spe- 
cific heat and entropy, for some specific cases. In partic- 
ular, we focus on the most important (n) ~ 1 situation. 
Our data can be used in two ways: (i) to understand lim- 
its of applicability of the semi-analytic approach (with 
calculated effective parameters) discussed above, and (ii) 
to have reference first-principle curves for more refined 
numerical analysis. Unfortunately, a direct simulation of 
a realistic case in the trap, i.e. with similar number of 
particles as in experiments, is still a challenging problem 
though simulations of about 10 5 particles or more at low 
temperature seem feasible in near future. 

The results are organized as follows: in Subsection 



A we probe the limits of applicability of semi-analytic 
predictions in the Mott state. In Subsection B we cal- 
culate the entropy of the Bose-Hubbard model, compare 
it to the initial entropy (i.e. before the optical lattice 
is turned on) and estimate the final temperature. We 
consider both a homogeneous system in the MI and SF 
states and a system of N ~ 3 • 10 4 particles in a trap. 



A. Comparison with low T semi-analytic 
predictions 

Away from the tip of the lobe, in the MI state, semi- 
analytic predictions are reliable provided the tempera- 
ture is low enough, i.e. r< A. In other words, there 
exist a range of temperatures , defined as T < const ■ A, 
where the quasi-particle excitations can be successfully 
described as a non-interacting classical gas (see Eq. (0). 
The value of the constant depends on A, as the two fol- 
lowing examples demonstrate. 

Let us first consider larger gaps, e.g. A ~ 200i, for 
which we have found that the low temperature analytic 
predictions reproduce numerical data very well. In Fig. 
|7| we plot the energy per particle in the low temperature 
regime for A = 181.64. The analytic prediction from Eq. 
([7]), where, at any given temperature, the chemical po- 
tential has been chosen by setting equal the total number 
of particle and hole excitations (as it is done for intrinsic 
semiconductors), is reliable up to temperatures T < 35t, 
i.e. T < 0.15A. In the inset we plot the average num- 
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FIG. 8: (Color online). Energy (left), specific heat (center) and entropy (right) per particle at t/U = 0.025 (MI ground state) 
and unity filling factor. On the left, solid circles refer to prime data (error bars within symbol size), the data were taken for 
linear system sizes L=10 and L—20. Within the error bars, we are not able to resolve any finite size effect. Solid lines in all plots 
are obtained from spline-interpolated data for energy, with subsequent analytic differentiation/integration of the interpolation 



ber of particle-hole excitations. This number increases 
rapidly with temperature justifying the grand-canonical 
calculation for the quasi-particle gas (at fixed total num- 
ber of particles). The quasi-particle number density is 
~ 5% at T ~ 35t. Apparently, for higher temperatures 
the ideal gas picture is no longer valid as it crosses over 
to that of the strongly correlated normal liquid. We con- 
clude that for large enough gaps and T < 0.15A, one 
can rely on low temperature analytical predictions to do 
thermometry. 

For smaller gaps, instead, we do not find any inter- 
esting region (i.e. where temperature effects are visible), 
for which the classical description is valid. In Fig. [8] we 
show the energy per particle as a function of temperature 
for t/U = 0.025 (the groundstate is MI with the energy 
gap A = 18.35t). To get the specific heat and entropy, 
we first use spline interpolation of the energy data points 
to obtain a smooth curve E(T). The specific heat is then 
obtained by differentiating the spline. The maximum in 
the specific heat is reached when temperature is about 
half the energy gap. The entropy has been calculated 
by numerical integration of cy/T. In order to see any 
temperature effect one has to go as high as T ~ 2.5t; at 
these temperatures the classical description is already no 
longer applicable and one has to rely on numerical data 
to do thermometry. 



B. Loading the optical lattice: estimate of T^ Rn > 
from entropy matching 



The standard approach to convert results obtained 
for a homogeneous system into predictions for systems in 
external fields is the so-called local density approxima- 
tion (LDA), which is actually a local chemical potential 
approximation when the density at the site i is identi- 
fied with the density of the homogeneous system with 



the chemical potential 



(off) 

M = M 



V{i) 



(11) 



In strongly interacting regimes with a short healing and 
correlation length, the LDA approach can be easily jus- 
tified in most cases (critical regions of phase transitions 
excluded). In Ref. [2l[, the authors directly compare 
simulation results for ID and 2D harmonically trapped 
systems with LDA predictions based on known homoge- 
neous system phase diagram. As expected, the density 
profiles differ only at the MI-SF interface, and we believe 
that the same will be true for the 3D case which is more 
"mean-field-like" . 

When the semi-analytic predictions are reliable (see 
Fig. [7J), one can use numerical results for the effective 
masses and gaps to calculate the entropy of the homoge- 
neous quasi-particle gas with the tight-binding dispersion 
relation. The entropy is given by: 



S = - 



V 



(2tt) 3 



3 . 9[fi+(k) + fi_(k)] 



where 



n±(k) 



T In 1 



dT 



exp • 



e ±(k) 
T 



(12) 



(13) 



As an example, consider a uniform weakly interact- 
ing Bose gas (WIBG) of 87 Rb with the gas parameter 
nal ~ 10 -6 , which is loaded into an optical lattice with 
A = 840nm and t/U — 0.005. At low enough temper- 
ature, T < 0.3T C , one can calculate the initial (prior to 
imposing the lattice) entropy of the system using the Bo- 
goliubov spectrum. In Fig. [5] we plot the entropy per unit 
volume before and after the optical lattice is imposed. 
The bottom x axis is temperature in units of the critical 
temperature of the WIBG, the top x axis is temperature 
in units of t, the hopping matrix clement. The dashed 
line is the entropy of the WIBG, the solid and dashed- 
dotted lines represent the entropy of the Bose-Hubbard 
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FIG. 9: (Color online). Entropy per particle at t/U = 0.005, 
unity filling factor and linear system size L = 20. The dashed 
line is the entropy of the uniform WIBG. The solid line is 
the result of analytical derivation/integration of numerical 
data for energy and the dashed-dotted line (barely visible) 
is the analytical prediction, Eq. (I12p . where, at each temper- 
ature, the chemical potential has been fixed by the condition 
of having equal number of particle and hole excitations. If the 
system was initially cooled down to T m = 0.25T C , the final 
temperature is T fin = 22t and nearly a hundred of thermally 
activated particle-hole excitations are present in the final state 
(see inset in Fig. [TJl . 

model calculated starting from the numerical results of 
Fig. Eland analytical predictions of Eq. (fT2"f . respectively. 
If the system was prepared at T ~ 0.25T C , the final tem- 
perature would be T ~ 22t. Fig. shows that at this 
temperature the system is quite far from its ground state 
and the number density of thermally activated particle- 
hole excitations is ~ 1%. The circumstances of this kind 
become crucial if one is to use the system in quantum in- 
formation processing. This example is also illustrative of 
how numerical data can be used to suggest appropriate 
initial conditions. 

Now we turn to a more realistic case of confined 
system and use LDA to convert results for the uniform 
system into predictions for the inhomogeneous one. In 
what follows, we consider a gas of N ~ 3 • 10 4 87 Rb 
atoms, magnetically trapped in isotropic harmonic po- 
tential of frequency 27r60 Hz. Experiments with such 
number of particles were recently performed [9j|. With 
this geometry, the parameter rj — 1.57(A rl / 6 a s /a/ lo ) 2 / 5 
(see Ref. (22J) is ~ 0.33, which is a typical value in current 
experiments. For temperatures in the range fi < T < T c , 
where T c is the critical temperature of the harmonically 
trapped ideal gas, one can accurately calculate energy 
using the Hartree-Fock [23| mean field approach [24[: 

^ = ffiWkl-i) 2/5 (5 + lM 3 ), (14) 



where t is the reduced temperature T/T c . At very low T 
(T < /i), Eq. (P3J misses the contribution coming from 
collective excitations. We are interested in initial tem- 
peratures T ~ 0.2 — 0.3T C , which are feasible in current 
experiments (25[. Starting from Eq. (fT4")l . we calculate 
the entropy of the BEC initially prepared in the magnetic 
trap. After the optical lattice is adiabatically turned on, 
the magnetic potential provides a scan over the chemical 
potential of the homogeneous system [see Eq. 

A direct comparison with experiments at fixed num- 
ber of particles would require to calculate n(T) from the 
normalization condition. At low temperatures, one ex- 
pects the dependence of the chemical potential on tem- 
perature to be weak (this will be confirmed by direct 
simulations, see below). For simplicity, we fix /i at a 
value corresponding to N = 30 3 trapped atoms in the 
first Mott lobe, at zero temperature. From this point we 
proceed in two directions. On one hand, we analytically 
calculate the low temperature contribution to energy and 
entropy arising from particle and hole excitations in the 
trapped MI state. On the other hand, we directly sim- 
ulate the thermodynamics of the inhomogeneous system 
at a fixed chemical potential. The results are shown in 
Fig. 1 101 where we plot the energy per particle, counted 
from the ground state. The solid circles are data from 
the simulations (error bars are plotted), the solid line is 
the (analytically calculated) contribution of the particle 
and hole excitations. The inset shows the low tempera- 
ture region. A large mismatch between the two results 
indicates that the main contribution to energy is given 
by the liquid at the perimeter of the trap. At zero tem- 
perature, there are about N ~ 29000 particles in the 
trap, 7% of which are not in the MI state (recall that 
fj, has been determined by placing 30 3 particles in the 
MI state). Simulation results show that, in the range 
of temperatures considered, the total number of parti- 
cles increases by 0.7% at most, which confirms the weak 
temperature dependence of the chemical potential. In 
addition, for T = 8t, we performed a simulation with N 
fixed at the groundstate value. The energies per particle 
in the canonical and grand canonical simulations differ 
by 0.3% only, therefore we proceed calculating the en- 
tropy in the canonical ensemble and compare it with the 
initial entropy of the system, at a fixed particle number. 

We are in a position to address the question of what 
is the final temperature of the system after the optical 
lattice is turned on and the final state is MI with the 
exception of a small shell at the trap perimeter. In Fig.fTTI 
we plot the entropy of the trapped WIBG with (solid 
line) and without (dashed line) the optical potential. If 
the initial system is cooled down to temperatures T( m ) ~ 
0.25T C , see, e.g., Ref. [25[, the final temperature will be 
T( fin ) = (2.35 ± 0.30)i.' 

With the initial conditions considered in this exam- 
ple, what is the final state of the liquid at the perimeter 
of the trap? Before answering this question, we would 
like to recall that, along the MI-SF transition lines, the 
critical temperature for the normal-to-superfluid transi- 
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FIG. 10: (Color online). Energy per particle at t/U = 0.005, 
fi — 116. 5t and trap frequency ui = 2n60Hz. Solid circles 
are numerical data, the solid line is the energy of particle 
and hole excitations in MI deduced from Eq. J7J|. The inset 
shows a zoom of the low temperature range. At T ~ 2.3t the 
contribution of MI excitations to energy is ~ 10% only. 



tion is zero. The transition temperature increases as one 
moves away from the border of the Mott lobes (lower- 
ing fi at fixed t/U in our case) and the quasi-particle 
density increases until it reaches its maximum at about 
n m 1/2 and then decreases. The maximum T c can 

be estimated from the ideal Bose gas relation = 
3.313n 2 / 3 /m* = 4.174t, with n = 0.5 and m* = l/2i, 
but interaction effects are likely to reduce this value. 
We have performed simulations at half filling factor and 
fixed t/U — 0.005, and found the critical temperature to 
be T c (n = 1/2) = 2.09(l)i. As a consequence, for the 
chosen initial conditions, we can conclude that the final 
state of the liquid at the perimeter of the trap is normal 
and it gives the main contribution to the entropy. For 
such low final temperature, the contribution to the en- 
tropy per particle due to thermally activated excitations 
in the MI state is only 10%. The largest chemical po- 
tential is, in fact, deep in the first Mott lobe, and the 
energy required to introduce an extra particle or hole is 
much larger than T. Most excitations are located at the 
perimeter of the Mott state in a narrow shell of radius R 
and width - 0.05i?. 

Retrieving the same information for experiments us- 
ing a larger number of particles, e.g. 10 5 10 6 , by di- 
rect simulation is still computationally challenging. In 
order to use LDA, one should study the uniform sys- 
tem, scanning through the chemical potential. As our 
last example, we consider a uniform system which is in 
the correlated SF ground state. Fig. [T^] shows data for 
E, cy, and S for t/U = 0.0385 and unity filling factor, 
close to the MI-SF transition. The system stays in its 
ground state for T <C 2t; in finite systems, the energy 
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FIG. 11: (Color online). Entropy per particle for the same 
system as in Fig llOl The dashed line is the entropy of the 
trapped WIBG, before the optical lattice is loaded; the solid 
line is the final entropy. If the system was initially cooled 
down to T (in) = 0.25T C , the final temperature is T (fln) = 
(2.35 ± 0.3)t and the liquid at the perimeter of the trap is 
normal (see text). 

of the lowest mode is finite: E m i n = cp m in, with c « 6t 
and p m i n = 2tt/L). The specific heat and entropy are 
calculated as described for Fig. [8] We were not able to 
resolve the SF-NL transition temperature from this set of 
data alone: numerical data corresponding to system sizes 
L = 10 and L — 20 overlap within error bars and we did 
not see any feature at the critical temperature. This is 
not surprising since the specific heat critical exponent a 
is very small, and it is thus very difficult to resolve the 
singular contribution and finite size effects in energy and 
specific heat. For a system of linear size L, the singular 
part of the specific heat, Cy , can be written as 

4\t,L) = C /v fcWL) = L a ^g c {tL^) , (15) 

where t = {T-T c )/T c , a w -0.01, v = (2-a)/3, £ is the 
correlation length and f c {x) and g c {x) are the universal 
scaling functions. At the critical point, finite size effects 
for the two system sizes considered are ~ 1%, within 
error bars. 

The critical temperature was extracted from data 
for the superfluid stiffness. The scaling of the superfluid 
stiffness at the critical temperature is n s oc \t\". This 
allows one to accurately estimate the critical temperature 
from 

n s (t,L) = C'm/L) = L- 1 9s (tL 1/lJ ) , (16) 

by plotting n s L vs. T as shown in Fig. 1131 From the 
data taken for system sizes L = 5, L = 10, L — 20, we 
estimate the critical temperature to be T c = 3.25(l)t. 
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FIG. 12: (Color online). Energy (left), specific heat (center) and entropy (right) per particle at t/U = 0.0385 (SF ground state) 
and unity filling factor. On the left, solid circles refer to prime data (error bars within symbol size). Data were taken for system 
size L=10 and L=20. Within error bars, we are not able to resolve any finite size effect. Solid lines in all plots are obtained 
from spline-interpolated data for energy, with subsequent analytic differentiation/integration of the interpolation curve. 




FIG. 13: (Color online). Finite size scaling of the super- 
fluid stiffness. n s L vs. T/t for system size L=5 (circles), 
L=10 (squares), L=20 (triangles). We estimate the critical 
temperature to be T c — 3.25(l)t, already nearly half the non- 
interacting gas value. 

At this temperature, the entropy per unit particle is 
~ 0.195, or, translating to entropy density in physical 
units, 3.6T0~ 5 JK -1 m~ 3 , which corresponds to an initial 
temperature ~ 0.35Tc \ Therefore it seems plausible to 
reach T c experimentally. 

IV. CONCLUSIONS 

We have performed quantum Monte Carlo sim- 
ulations of the three dimensional homogeneous Bose- 
Hubbard model. We were able to establish the phase 
diagram of the MI-SF transition with the record accu- 
racy ~ 0.1% and determine the size of the fluctuation 



region in the vicinity of the diagram tip where universal 
properties of the relativistic effective theory can be seen. 
Comparison with the strong-coupling expansion shows 
that the latter works well only for sufficiently large insu- 
lating gaps A > 6i outside of the fluctuation region. Wc 
have studied the effective masses of particle and hole ex- 
citations along the MI-SF boundary. Our results directly 
demonstrate the emergence of the particle-hole symme- 
try at the diagram tip, and provide base for accurate 
theoretical estimates of the MI thermodynamics at low 
and intermediate temperatures. 

We have studied thermodynamic properties of the 
superfluid and insulating phases at fixed particles num- 
ber for the uniform case. These data can be used to make 
predictions for the inhomogeneous system using the lo- 
cal density approximation. We have shown that for large 
enough gaps the low temperature analytical predictions 
agree with numerical data. By entropy matching, we 
have calculated the final temperature of the system (af- 
ter the optical lattice is adiabatically loaded), in the uni- 
form and magnetically trapped system, at t/U = 0.005. 
We have performed direct simulations of a trapped sys- 
tem, using typical experimental values for the magnetic 
potential and number of particles. For the initial con- 
ditions considered, we found the final temperature and 
demonstrated that the main contribution to the entropy 
comes from the liquid at the perimeter of the trap. Wc 
have calculated the normal-to-superfluid transition tem- 
perature at the half filing and concluded that the liquid 
at the perimeter is normal. 
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